Center for Turbulence Research 
Annual Research Briefs 2003 


173 


New type of the interface evolution in the 
Richtmyer-Meshkov instability 

By S. I. Abarzhi and M. Herrmann 


1. Motivation and objectives 

When a shock wave passes an interface between two fluids with different values of the 
acoustic impedance, the misalignment of the pressure and density gradients results in 
a growth of the interface perturbations and causes the development of the Richtmyer- 
Meshkov instability (RMI) (Richtmyer 1960; Meshkov 1969). The instability produces 
with time the turbulent mixing of the fluids, which controls many physical and techno- 
logical processes, such as inertial confinement fusion, supernova explosion, and impact 
dynamics of liquids (Kull 1991). Reliable description of the turbulent mixing is the basic 
objective of studies of RMI (Kull 1991; Schneider et al. 1998). 

Observations report the following evolution of the Richtmyer-Meshkov instability. In 
the linear regime, the light fluid accelerates impulsively the heavy fluid, and a small 
amplitude perturbation of the fluid interface grows linearly with time (Richtmyer 1960; 
Meshkov 1969). The acceleration value is determined by the shock-interface interaction, 
which is essentially non-local and results in a baroclinic production of vorticity (Haan 
1991; Velikovich & Dimonte 1996; Vandenboomgaerde et al. 1998; Wouchuk 2001). In 
the non-linear regime, a structure of bubbles and spikes appear (Pavlenko et al. 2000; 
Chebotareva et al. 1999; Jacobs & Sheeley 1996; Bonazza & Sturtevant 1996). The light 
(heavy) fluid penetrates the heavy (light) fluid in bubbles (spikes). The bubbles decelerate 
and the spikes move steadily. Small-scale structures appear on the side of evolving spikes 
due to shear, and the direct cascade of the fluid energy may occur (Matsuoka et al. 2003). 
For a finite-amplitude perturbation, the fluid energy may be transferred also to larger 
scales (Alon et al. 1995; Oron et al. 2001). Eventually, a mixing zone develops. In the 
chaotic regime, the bubbles and spikes decelerate, and their positions are described by 
power-law time-dependencies with exponents determined by the density ratio (Schneider 
et al. 1998; Dimonte 2000). 

The dynamics of the Richtmyer-Meshkov instability is far from being completely un- 
derstood. For a long time, the theoretical models and numerical simulations failed to 
predict the growth-rate in the linear RMI observed in experiments (Holmes et al. 1999). 
Only recently, a rigorous theory has accounted for the non-local character of the inter- 
face dynamics in the case of strong and weak shocks and provided therefore for impulsive 
models the correct value of the acceleration (Wouchuk 2001). To describe the nonlin- 
ear RMI, several models have applied a single-mode approximation, which presumed 
locality of the bubble (spike) dynamics (Alon et al. 1995; Goncharov 2002). Some other 
models have used an empiric equation with adjustable parameters to balance “drag” 
and “inertia” in the flow (Oron et al. 2001; Dimonte 2000). All these models however 
cannot explain the observations and remain subjects for controversy (Dimonte 2000). 
There is a strong need in a formal theoretical approach and in experiments and simula- 
tions with systematic variation of parameters and improved diagnostics. In this work we 
suggest analytical and numerical solutions describing the nonlinear coherent dynamics 
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Figure 1 . Interface definition. 


of two-dimensional Richtmyer-Meshkov instability for fluids with a finite density ratio. 
Our results predict new properties of the nonlinear evolution of RMI, explain existing 
experiments, and identify new sensitive diagnostic parameters. 

2. Governing equations 

The dynamics of the Richtmyer-Meshkov instability are governed by a system of con- 
servation laws, which are compressible two-dimensional Navier-Stokes with the initial 
conditions and the boundary conditions at the fluid interface, 


| + v. W = o 

(2.1) 

dpv — 

-gj- + V • (pvv) + 'S/p + V • t = 0 

(2.2) 

q +S7 ■ ([pE + p]v) + S7 ■ q + t : S7v = 0. 

(2.3) 

Here, p denotes the density, v the velocity vector, p the pressure, r the stress tensor, E 
the total energy, and q the heat flux vector. All quantities are either in the heavy gas, 
denoted by the subscript h in the following, or in the light gas, denoted by the subscript 
l. The above system of equations is closed by the ideal gas law, 

p = pA>T , 

(2.4) 

with <I> the gas constant and T the temperature. 

The interface T located at z*(x,t) separating heavy from light gas, 
level set scalar G, see Fig. 1. Defining 

is described by a 

G(x , t) = Go = const , 

(2.5) 


with G(x,t) < Go in the light gas and G(x,t) > Go in the heavy gas, compare Fig. 1, 
an evolution equation for the scalar G can be derived by simply differentiating Eq. (2.5) 
with respect to time, 

r)G 

— + wVG=0. (2.6) 

at 

This equation is called the level set equation (Osher & Sethian 1988). It is easy to see 
that Eq. (2.6) is independent of the choice of G away from the interface. However, to 
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facilitate the numerical solution of Eq. (2.6), G is chosen to be a distance function away 
from the interface, 


|VG| 


G^G 0 


= 1 . 


(2.7) 


Using the level set scalar, geometrical properties of the interface, like its normal vector 
n or curvature Q, can be easily calculated, 


VG 


k = V • n . 


( 2 . 8 ) 


There is no mass flow across the moving interface, and the normal component of 
velocity, pressure and temperature are continuous at the fluid interface: 



(2.9) 

Vt,h | p = V t,l | p 

(2.10) 

Ph r =Pi r 

(2.11) 

Tl 

II 

(2.12) 


The flow has no mass sources, and the boundary conditions at the infinity close the 
set of the governing equations 


V *\z=+oo = V ‘\z=-oo =0 - ( 2 - 13 ) 

Initially, the fluid interface is slightly disturbed with a small amplitude co-sinusoidal 
perturbation, z*(x,t = 0) ~ (1/fe) cos(fcr), where k = 2 tt/\ is the wave-vector, and A 
is the spatial period in the x-direction. The perturbation should be symmetric, in order 
a stable coherent structure of the bubbles and spikes to occur. The density ratio or the 
Atwood number A = (p/, — pi)/(ph + Pi) is a determining factor of RMI dynamics. 


3. Non-local theoretical solutions 

Based on the observations, we divide the fluid interface into active regions (small scales) 
with intensive vorticity, and passive regions (large scales) which are simply advected. If 
the energy cascades are not extensive (i.e., the fluid densities are not very similar, the 
perturbation amplitude is small, and the initial shock is weak), then, a considerable part 
of the fluid energy concentrates in the large-scale coherent motion with V • v h ^ = 0 
and V x v^i) = 0. To find the nonlinear solutions, describing the dynamics of the 
bubble front in a vicinity of its tip, we reduce the Eqs. (2.1), (2.2), (2.9), and (2.11) to a 
local dynamical system. All calculations are performed in the frame of reference moving 
with velocity v(t) in the 2 -direction, where v(t) is the velocity at the bubble tip in the 
laboratory frame of reference. For the large-scale motion v^i) = and we expand 

the potential $h(i) as a Fourier series, 

OO 

<!>/, = ^2 $ m (t) (cos(mfcx) exp(-mkz)/mk + z) 

m = 1 
oo 

d >; = ^2 < Fm(i)(cos(mfcx) exp(mkz)/m — 2 ) 

m—1 

For x « 0 the interface has the form 2 * = (N{t)x 2N , where £i(t) < 1 is the 


(3.1) 

(3.2) 
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principal curvature at the bubble tip, and N is the order of approximation. Substituting 
these expressions in the governing equations, taking the first integral of Eq. (2.2), and 
re-expanding Eqs. (2.9) and (2.11) for x ss 0, we derive a system of ordinary differential 
equations for the surface variables Civ(t), and the moments M n (t) = 5Zm=i ®m{t){km) n 
and M n (t) = $ m {t){km) n where n is integer. For N = 1, one has from Eqs. (2.9)- 

(2.13) respectively 


Cr = 3CiMi + M 2 / 2 = 3CiMi - M 2 /2, (3.3) 

(Mr/2 + Ci M 0 - M\l 2 - Ci)p h = (Mr/2 - Ci^o - M?/2 - Ci)«> (3-4) 

M 0 (t) = -M 0 (t) = -v(t) . (3.5) 


The system (3.3)-(3.5) describes the local dynamics of the bubble as long as the period 
of the structure is invariable, and the energy cascades are not extensive. The presentation 
in terms of moments allows one to account for the effect of higher-order correlations. The 
time-scale in Eqs. (3.3)-(3.5) is r = 1/kvo, where vq is the absolute value of the initial 
velocity. For t/r <C 1, the dynamical system (3.3)-(3.5) has regular asymptotic solutions 
with time independent surface variables Civ and with velocity v and moments M n , M n 
decaying as 1/t. 

It is easy to see that the local system (3.3)-(3.5) cannot be satisfied in a single-mode 
approximation. One may derive, for example, a single-mode solution of the Layzer- 
type, which conserve mass, momentum and has no mass sources. For this bubble Cl = 
Cl = —Ak/6 and velocity v = vl = (1 — A 2 /3)/Ak, however the solution does not 
satisfy Eqs. (2.9) and (3.3) and requires mass flux across the interface. To reproduce 
the parameters of the drag model (Oron et al. 2001) with Ci = Cd = — fc/6 and 
v = vd = (1 + A/3)/k(l + A ), one should violate Eqs. (2.13) and (3.5) and intro- 
duce a source of mass of the light fluid (Goncharov 2002). Obviously, these solutions are 
unphysical, because they violate the conservation laws. 

To obtain regular asymptotic solutions describing the nonlinear dynamics of the bubble 
front, one should account for the non-local properties of the flow that has singularities. 
The singularities determine the interplay of harmonics in the global flow and the local 
dynamical system, and affect therefore the shape of the regular bubble. We find a con- 
tinuous family of regular asymptotic solutions for Eqs. (3.3)-(3.5), parameterized by the 
principal curvature at the bubble tip, and choose the fastest solution in the family as the 
physically significant one. For TV = 1 the bubble velocity as the function of the bubble 
curvature is 


v = 


3 (1 + Afa/k) - 12A(Ci A) 3 ) 
2 kt (A - 4(Ci/fc) + 4A(Ci/fc) 2 ) 


(l-4(CrA) 2 ), 


(3.6) 


where ( cr < ( < 0. For the family solutions (3.6), the interplay of harmonics is well 
captured, the higher order corrections for the velocity and lowest-order amplitudes are 
small, the solutions converge, yet, most of them are unstable. The fastest stable solution 
in the family (3.6) is the solution with 


Ci = (a, = 0 v = va = 3/2 Akt 


(3.7) 


The foregoing theoretical results suggest the following evolution of the bubble front in 
the Richtmyer-Meshkov instability. In the linear regime of RMI, the bubble curvature and 
velocity change as ~ t ; in the weakly non-linear regime, the curvature reaches an extreme 
value, dependent on the initial conditions and the Atwood number; asymptotically, the 



New type of the interface evolution in the Richtmyer- Meshkov instability 177 

bubble flattens and decelerates. For A < 1 the bubbles move faster than those for A = 1, 
and for all A the bubbles flatten asymptotically. The flattening of the bubble front is 
a distinct feature of RMI universal for all A. It follows from the fact that RM bubbles 
decelerate. As A ~ 0, the velocity v ~ oo and this suggests that for fluids with similar 
densities the bubble velocity has a much faster time-dependence, such as t a where —1 < 
a < 0, in a qualitative agreement with experiments of Jacobs & Sheeley (1996), where 
for A ~ 0 the dependence 1 /t was shown to underestimate the velocity data. 

4. Numerical method 

In our numerical simulations, the Navier-Stokes equations (2.1)-(2.3) are solved using 
a hybrid capturing-tracking scheme, originally proposed by Smiljanovski et al. (1997) for 
deflagration waves. When applying this scheme to the RMI, its key idea is to explicitly 
track the location and motion of the interface between the light and the heavy gas by 
the level set equation (2.6), whereas all other fluid phenomena like shocks and expansion 
fans are captured. The main advantage of this approach is that while the simplicity 
and robustness of standard capturing schemes can be retained, the interfacial processes 
are described with accuracy comparable to standard tracking schemes. In the following 
two sections, the hybrid capturing-tracking scheme is summarized briefly. The interested 
reader is referred to Smiljanovski et al. (1997) for a more detailed description. 

4.1. In-cell-reconstruction 

In finite volume schemes, the cell value of a conserved quantity U 1 ’ 3 is defined as the 
volume average of that quantity, 

UiJ = VU i vij U(x') d x\ ( 4 . 1 ) 

averaged over the cell volume V 1 ' 3 . Assuming piecewise constant distributions of U for 
each fluid within each cell, Eq. (4.1) reduces to 

U 1 ' 3 = aU^ 3 + (1 - a)U\ 3 , (4.2) 

where a is the heavy gas cell volume fraction that can easily be calculated from the level 
set scalar G, 

a=^. j^H{G{x')-G 0 )dx', (4.3) 

with H the Heavyside function. 

The key idea of the in-cell-reconstruction scheme is to reconstruct both U ff and U\'° 
from U 1 ' 3 , see Fig. 2, and to use only the reconstructed two states to calculate cell face 
fluxes in cells containing part of the interface. Combining the jump conditions of U across 
the interface, Eqs. (2.9) - (2.12), with (4.2) and Eqs. (4.3), the cell values U]’ 3 and TJ\ 3 
in each cell containing part of the interface can be reconstructed from U 1 ’ 3 . 

4.2. Cell update 

The numerical solution is evolved in time using an operator splitting technique (Strang 
1967), 

U n+1 = C* At/2 C* At/2 V At C z At/2 C* At/2 U n , (4.4) 

where C denotes the convection operator, T> the diffusion operator, U n the cell volume 
averaged solution at time t n , and U n+1 the cell volume averaged solution at t n+1 = 




a 

Figure 2. In-cell-reconstruction, 
time level t" flux upda te 



Figure 3. Flux update in cells containing part of the interface. 


t n + At. In order to calculate the correct update, all gradients in the individual operator 
steps must be calculated using only cell values of the same fluid type. 

In the diffusion operator T>, we calculate all gradients using 2nd order central differ- 
ences. The corresponding stencils across the interface thus involve the reconstructed cell 
values Uh respectively Ui plus one additional adjacent cell value on the opposite side of 
the interface. This cell has to be first transformed into the corresponding matching state 
using Eqs. (2.9)-(2.12) before the gradients are evaluated. In this operation, our method 
resembles the ghost fluid approach (Fedkiw et al. 1999). 

Since the convection operator C is split into each spatial direction, we will only focus on 
the ^-direction operator in the following. If a cell (i,j) and its directly adjacent neighbors 
do not contain part of the interface, U n ' l ’ J is advanced in time by 

jjn+l,i,j _ jjn,i,j _|_ ^pi-l/2,j _ pi+l/2,j^j . ( 4 . 5 ) 

In our numerical method we solve cell face Riemann problems using a 2nd order wave 
distribution algorithm due to LeVeque (1990) and formally recast the individual wave 
contributions in the form of cell face fluxes _F* _1 / 2j and _F ,J+1 / 2 >T However, special care 
must be taken, if (i, j) contains part of the interface at t n or t n+1 , as shown in Fig. 3. 

To ensure the correct flux calculation, individual cell face fluxes for both the heavy and 
the light gas must be calculated using only (reconstructed) quantities of the respective 
fluid, 

77 U— 1/2, i Tpi— 1/2, J f rri— 1/2, J j-j-i— 1/2, j 

* h ~ * h \ h ~ ’ ^ h+ 

pi— 1/2, j pi— 1/2, j ^jji— 1/2, j jji— 1/2, j 


(4.6) 

(4.7) 
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Here, the subscripts — and + denote the left-hand side respectively right-hand side 
approximation of U‘ h p 2 ’- 7 . Then, the average cell face flux F l ~ 1 ^ 2 ’ 3 is 

F *- 1/2 0 = jji-l/2,j pi-l/2,j + (j _ pi-1/2, j^pi-1/2,3 , ( 4 . 8 ) 

with f3 the heavy gas cell face fraction, see Fig. 3. Note that because our convection 
operator is based upon a wave distribution scheme, F 1-1 ^ 2 ’ 3 only takes those wave con- 
tributions into account that arise from the cell face Riemann problem. Since the interface 
itself constitutes a separate wave, its contribution to the change of U 1 ’ 3 has to be taken 
into account additionally, 

AC/*’- 7 = Aa c ( U\ 3 - U\ 3 ) + 

Aa L (t/p 1 ’- 7 - C/p 1 ’- 7 ) + Aa R ([/p 1 ’- 7 - C/p 1 ’- 7 ') . (4.9) 

Here, Aa c is the change of the heavy gas cell volume fraction due to the movement of 
the interface within the cell (i,j) itself, and Aa L and Aa R are the contributions due to 
movement from adjacent cells to the left respectively right, see Fig. 3. Note that Eq. (4.9) 
implies that the global conservation property of the numerical scheme is now dependent 
on the accuracy with which the individual Ao are calculated and hence relies on the 
accuracy of the level set method. 

Finally, combining Eq. (4.8) and Eq. (4.9) yields the update for U 1 ' 3 in cells that 
contain part of the interface 



AcP (r/p - i/p) + An'- (ut 1 ’ 3 - ur 1 ’ 3 ) + 

Aa R (c/p 1 ’- 7 — t/p 1 ’- 7 )} (4.10) 

One of the major advantages of using Eq. (4.10) is the fact that the CFL-condition can 
be based on the grid size Aa; of the underlying grid, because cell updates are performed 
only on the cell volume averaged quantities. 


5. Results 

We present the numerical results of the RMI calculated for a shock Mach number of 
Ma =1.2 and two different Atwood and Reynolds numbers. The fist case corresponds 
to the SFg/Air experiment of Benjamin et al. (1993) with A = 0.67 and Re = 11442, 
whereas in the second case A = 0.9 and Re = 6977. The Reynolds numbers are calculated 
using the bubble velocity V in the laboratory frame of reference and the wave length A. 
The initial interface of wave length A = 0.0375 m and amplitude ao = 0.0024 m is located 
at z = 0 m. All simulations are performed in a [—1.1m, 0.3 m] x [—0.01875 m, 0.01875 m] 
box resolved by 2352 x 63 equidistant cartesian grid cells. 

Figure 4 depicts the temporal evolution of the interface shape for A = 0.67, whereas 
Fig. 5 shows the interface shapes for A = 0.9. In both cases, a spike of heavy gas is 
formed that penetrates ever further into the light gas. In the A = 0.67 case, the spike 
takes on the typical mushroom-like shape, whereas in the A = 0.9 the spike simply bulges 
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Figure 4. Temporal evolution of the interface for A=0.67. 
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Figure 5. Temporal evolution of the interface for A=0.9. 

A =0.67 A =0.90 




Figure 6. Shape of the bubble front for A=0.67 (left) and A=0.9 (right) plotted every 
A t/r = 10. Each interface is shifted by A a = —0.05. 


at the end. This difference is most likely due to the employed numerical resolution that 
is not sufficient to resolve the apparently smaller scale spike mushroom structure in the 
A = 0.9 case. 

Figure 6 shows the temporal evolution of the bubble front for both cases. The flattening 
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Figure 7. Bubble velocity V in the laboratory frame of reference for A = 0.67 (left) and 

A = 0.9 (right). 




Figure 8. Bubble curvature ( for A = 0.67 (left) and A = 0.9 (right). 


of the bubble front that is predicted by non-local theory is clearly visible in both cases, 
but is more pronounced in the A = 0.67 case. Also, the A = 0.67 case exhibits some tiny 
oscillations of the bubble front at later times that are not visible in the A = 0.9 case. 

The bubble velocity V in the laboratory frame of reference is depicted in Fig. 7. Initially 
being at rest, both bubbles are impulsively accelerated by the passing shock. As predicted 
by non-local theoretical analysis, the bubble velocity reaches a local maximum in the 
weakly non-linear regime and then asymptotically decelerates to a constant velocity. In 
the laboratory frame of reference the asymptotic bubble velocity is V = 66.8 m/s for 
A = 0.67 and V = 40.7 m/s for A = 0.9. 

Figure 8 shows the evolution of the bubble curvature (. Here, ( is calculated using a 
least squares fit of a circle of radius 1/|CI to all intersection points of the interface with 
the cell face within \x\/X < 0.12. As predicted by non-local theoretical analysis, in the 
linear regime, the bubble curvature changes linear with time. Then, in the weakly non- 
linear regime, the curvature reaches an extremum followed by an asymptotic flattening 
of the bubble, i.e. a decrease in curvature for both values of A. 

Finally, Figure 9 depicts the calculated bubble velocity v as a function of the absolute 
bubble curvature |£|. Initially, the bubble exhibits an abrupt acceleration that is due 
to the interaction with the passing shock, whereas the curvature remains roughly un- 
changed. In the linear regime, the bubble curvature increases with only gradual changes 
in the bubble velocity. This result is consistent with linear theory, Eqs. (3.3)-(3.5). In 
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Figure 9. Bubble velocity v as function of absolute bubble curvature |£| for A = 0.67 (left) 

and A = 0.9 (right). 


the weakly non-linear regime, the bubble velocity reaches a local maximum and then 
starts to decrease with almost constant curvature. Finally, the bubble curvature asymp- 
totically tends to zero while the bubble continues to decelerate, as predicted by non-local 
theory. However, depending on the case, a slightly different bubble behavior is evident. 
For A = 0.67, the bubble both flattens and decelerates simultaneously. For A = 0.9, 
the bubble deceleration occurs prior to any significant decrease of (, thus leading to a 
flattening of the bubble front with only gradual changes in the bubble velocity. 

In summary, these numerical results affirm the bubble behavior predicted by the non- 
local theory. 


6. Conclusion 

We performed systematic theoretical and numerical studies of the nonlinear large-scale 
coherent dynamics in the Richtmyer-Meshkov instability for fluids with contrast densi- 
ties. Our simulations modeled the interface dynamics for compressible and viscous fluids. 
For a two-fluid system we observed that in the nonlinear regime of the instability the 
bubble velocity decays and its surface flattens, and the flattening is accompanied by 
slight oscillations. We found the theoretical solution for the system of conservation laws, 
describing the principal influence of the density ratio on the motion of the nonlinear 
bubble. The solution has no adjustable parameters, and shows that the flattening of the 
bubble front is a distinct property universal for all values of the density ratio. This prop- 
erty follows from the fact that the RM bubbles decelerate. The theoretical and numerical 
results validate each other, describe the new type of the bubble front evolution in RMI, 
and identify the bubble curvature as important and sensitive diagnostic parameter. 
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